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The physical behaviour of a class of mesoscopic models for multiphase flows is analyzed in details 
near interfaces. In particular, an extended pseudo-potential method is developed, which permits to 
tune the equation of state and surface tension independently of each other. The spurious velocity 
contributions of this extended model are shown to vanish in the limit of high grid refinement and/or 
high order isotropy. Higher order schemes to implement self-consistent forcings are rigorously com- 
puted for 2d and 3d models. The extended scenario developed in this work clarifies the theoretical 
foundations of the Shan-Chen methodology for the lattice Boltzmann method and enhances its 
applicability and flexibility to the simulation of multiphase flows to density ratios up to 0(100). 

PACS numbers: 68.03.Cd,05.20.Dd,02.70.Ns,68.18. Jk 

INTRODUCTION 

The Lattice Boltzmann method [110, Q developed in the late 80's as an efficient and powerful way to simulate nearly 
incompressible hydrodynamics and its multiphase extensions 0, S S] represent one of the most successful mesoscopic 
techniques for numerical simulation of complex flows. 

Besides the mainstream application, namely complex macroscopic flows far from equilibrium, recent work is also 
hinting at the possibility that multiphase lattice Boltzmann methods may provide a new methodological framework 
for the description of fluid-solid interactions which play a crucial role for micro/nano-fluidic applications 0, Si- For 
example, the possibility to model slip boundary conditions and wetting properties [^, [l^, [h., \X, J^, [l^ has been 
recently achieved within the framework of the lattice Boltzmann equation. More detailed comparisons between the 



mesoscopic technique and atomistic molecular dynamics simulations |15l . Ilq have pointed out that lattice Boltzmann 
may become a method of choice for physical problems where supramolecular details play a major role. By supramolec- 
ular, we refer to situations which escape a purely continuum treatment, and yet, still exhibit sufficient universality to 
do away with a fully atomistic description. Arguably, a wide class of multiphase flows out of equilibrium falls within 
this class. 

All this looks promising, especially in view of recent experimental activity aimed at shedding some light on the rich 
and still largely unexplored territory of dynamical behaviour of liquids confined at (or below) millimetric scales. Impact 
of droplets on solid substrates, droplets breakup, capillarity instabilities and bouncing transitions, liquid fragmentation 
and water repellency on structured surfaces, are just but a few examples in point [ttI . [isl . ITol . |2a,[2l|,|23,[2|. 

Since the phenomenological description is not based on molecular details but only on average properties (for example 
surface tension, contact angle) mesoscalc modelling and numerical simulations would be extremely helpful to access 
time and space scales of direct experimental relevance. 

This is confirmed by recent numerical simulations for static behaviour [^, [l5| and also by some attempts to describe 
contact line motions ^24. 25, 26] and dynamical properties induced by heterogeneous wetting (ol. [lol|. 



These recent developments unquestionably set a compelling case for revisiting and extending some basic theoretical 
aspects of multiphase mesoscopic methods. In particular, the pseudo-potential approach introduced a decade ago 
by Shan and Chen (SC) 0, % to deal with non-ideal inhomogeneous fluids, represents one of the most successful 
outgrowths of the Lattice Boltzmann theory. It is worth noticing that non-ideal fluid behaviour can also be encoded 
a-priori by deriving lattice local equilibria directly from a free-energy functional 0]. This option leads to local 
equilibria with an explicit dependence on the density gradients, which cannot be readsorbed into a compact shift of 
the velocity, as it is the case for the pseudo-potential method 0, Q . The result is that the pseudo-potential method, 
albeit in-principle less rigorous, is very flexible and robust for practical and numerical purposes. 
Despite its undeniable success, this method has made the object of extensive criticism, the major objections being 
that surface tension is not tunable independently of the equation of state and that the interface dynamics is affected 
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by spurious currents near (curved) interfaces. 

In this paper, it is shown that both above hmitations can be hfted by moving to a mid-ranged pseudo-potential, i.e. 
by extending the spatial range of the pseudo-potential interaction. More specifically, it will be shown that (i) surface 
tension can be tuned independently of the equation state, by formulating a two-parameter version of the SC model 
with mid-range interactions, (ii) spurious currents near curved interfaces become vanishingly small in the limit of zero 
mesh-spacing and/or in the limit of an isotropic lattice. These developments help to put pseudo-potential methods 
a-la Shan-Chen on a solid theoretical basis. 



MEAN FIELD APPROACH: SHAN-CHEN MODEL AND ITS GENERALIZATIONS 



In this section, we briefly recall the main features of the lattice Boltzmann equation and the application to multiphase 
flow via the introduction of a pseudo-potential. The main goal here is to understand the corrections to the ideal-gas 
equations introduced by the presence of attractive pseudo-potential between Boltzmann kinetic populations. 

We start from the usual lattice Boltzmann equation with a single-time relaxation 27, 2^: 



fi{x + ciAt, t + At)- fi{x, t) = -— (fi{x, t) - fl"'\p, pu) ) + Fl (1) 



T 

where fi{x, t) is the kinetic probability density function associated with a mesoscopic velocity c;, t is a mean collision 
time (with At a time lapse), fi'^'^\p, pu) the equilibrium distribution, corresponding to the Maxwellian distribution 
in the continuum limit and Fi represents a general forcing term whose role will be discussed later in the framework of 
inter-molecular interactions. From the kinetic distributions we can define macroscopic density and momentum fields 
as [i,[2|: 



p{x)^Y.f^(^^ (2) 



pu{x) = ^cifi{x). 



(3) 



For technical details and numerical simulations we shall refer to the nine-speed, two-dimensional 2DQ9 model 0], 
often used due to its numerical robustness [29l |. The equilibrium distribution in the lattice Boltzmann equations is 
obtained via a low Mach number expansion of the continuum Maxwellian [s, 28J 
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(4) 



where 



= 1/3 and i = 1,2 = x,y runs over spatial dimensions. The weights W;'^''* are chosen such as to enforce 
isotropy up to fourth order tensor in the lattice 0, \2§^ . From the equilibrium distribution and the symmetry properties 
of c;, it immediately follows Q the kinetic second order tensor of the equilibrium distribution: 

I 

where, in the first term of the rhs, we recognize the well-known ideal-gas pressure tensor: 



Pi. 



(5) 



In order to study non-ideal effects we need to supplement the previous description with an interparticle forcing. This is 
done by choosing a suitable Fi in ([1]) . In the original SC model 0, Q , the bulk interparticle interaction is proportional 
to a free parameter (the ratio of potential to thermal energy), C/b, entering the equation for the momentum balance: 



F^ = -gbclYw{\ci\^)il){x,t)il}{x + ClAt.tYi 
I 



(6) 
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being w(|c;p) the static weights (w(l) = 1/3, w(2) = 1/12 for the standard case of 2DQ9 and ip{x^t) = ■>p{p{x,t) 
the (pseudo) potential function which describes the fluid-fluid interactions triggered by inhomogeneities of the density 
profile. The only functional form of the pseudopotential ijj{p) strictly compatible with thermodynamic consistency is 
V'(P) = p 0,[33- For purposes which will become apparent in the sequel, here we shall refer to the pseudopotential 
used in the original SC work namely 

^{x) = il-exp{-p{x))). (7) 

Note that this reduces to the correct form ip ^ p in the Hmit p <C 1, whereas at high density {p ^ 1), it shows a 
saturation. This latter is crucial to prevent density collapse of the high-density phases (note that the SC potential is 
purely attractive, so that a mechanism stabilizing the high-density phase is mandatory to prevent density collapse). 
In principle, other functional forms may be investigated, sometimes with impressive enhancement of the density 
ratios supported by the model js^l . 

In order to understand the corrections to the ideal-state equation ^ induced by the pseudo-potential, we need to 
define a consistent pressure tensor, P^, for the macroscopic variables: 

d,P,,^-F, + d,{c^p). (8) 

Upon Taylor expanding the forcing term and assuming hereafter At = 1, we obtain 

= -Gbijd^^ - y ^9,A^ + Oid^) (9) 

which is correctly translated into 

P^J - (^cip + Gt^^j'- + Qb^MA' + Gbji'A^j?j % - ^Gbcid.^d,^ + 0{d^). (10) 

Let us notice that there exists a sort of gauge-invariance in the definition of the pressure tensor, and pop is just 
one of these. In fact, while the term ipdiip is uniquely written as the gradient of V'^/2, the same is not true for the 
term ipdiAip. There are infinitely many tensorial structures that correspond to the same ipdiAip. On the other hand, 
from its very definition, it is clear that the tensor Pij is defined modulo any divergence-free tensor. However, it 
can be shown (see Appendix A) that all tensorial structures consistent with the forcing yield the same macroscopic 
surface tension and density profiles across the interface. Dispensing with consistency between the forcing term and 
the pressure gradient in the continuum, several choices for the pressure tensor can be proposed Hereafter, wc will 
stick to the requirement to have any of the possible gauge-invariant definition of the forcing and use the expression 
pO)) for all subsequent technical developments. 



In order to calculate the density profile for a flat interface in 2d whose dishomogeneities develop along a single 
coordinate, say y, we follow the mathematical details discussed in 0, [l3| and impose the mechanical equilibrium 
condition for the normal component Pyy of the above pressure tensor 



d P 







(11) 



with the boundary conditions that p(— oo) = pg and p(+oo) = pi , being Pg,pi the densities of the two phases. After 
some lengthly algebra (see [l^l for all details) one can show that the densities in the two phases are fixed by an integral 
constraint 
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where Pbuik 

r,2 



defines the equilibrium bulk pressure in one of the two phases, P; 



bulk 



cIpi 



c'Gb 



2 i'ipir = 

see also [3l|) because the 

latter is derived without imposing consistency with the forcing ([8]). Anyhow, for all practical purposes, the difference 
in the density ratios between the two versions is fairly negligible (see figure 7 of [lol|). 



c^Pg + ^^^V'(Pg)^- Let us notice that expression ([12]) is different from equations (25) of 



The surface tension also follows as the integral along a flat surface of the mismatch between the normal Pyy 
transversal Pxx component of the pressure tensor 32] : 



and 



(Pyy - Pxx)dy = - 



5b4 



\dyip\'^dy. 



(13) 
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Pseudo-potential with mid-range interactions 



It is immediately realized that, since in the SC model there is just a single free-parameter, Qb, it is impossible to 
tune density ratios (i.e. equation of state) and surface tension (i.e interface width) independently. In order to discuss 
this problem, let us go back to the expression of the forcing and consider possible generalizations thereof. The most 
immediate generalization of the standard SC model reads as follows: 

F, = -clY,wM)^i^)[SM^ + ci) + g2Mx + 2ci)]cl (14) 

where interactions up to next-nearest neighbors are explicitly enabled. The corresponding equilibrium pressure-tensor 
takes now the form: 

= (^clp + A,^^' + A,^\V4'\'+A2^^PA^py,,-^A2C%^l'^J^^ + Oi^'), (15) 

with Ai,A2 macroscopic constants related to Gi,G2 in (HI]): 

A2=gi+ 202 Ai=gi+ 802. (16) 

The surface tension now becomes 

/+00 A 4 r+oc 
[Pyy - P..)dy = Y / \^v^\^dy, (17) 

with the profile obtained applying the mechanical stability equation pT|) to fTS]) . Let us first notice that the above 
two-parameters couplings can be viewed as the first two terms of the expansion in terms of moments of the interacting 
potential: Ap = X]n=i° "■''^P' which is the lattice equivalent of the continuum virial expansion: Ap = J drr'PV{r), 
where V{r) is a general atomistic interaction potential. In principle one could also enlarge the spectrum of mid-range 
interactions but, for our purposes, it is enough to consider a two parameter coupling in (|14[) . Infact, in the case of 
equation (|15p we have an expression depending on the two free-parameters, Ai, A2, and on the functional shape of tp 
as a function of p. This opens up new degrees of freedom with respect to the standard SC formulation. First, let us 
fix the pseudo-potential shape to: 

V'(P) = - exp{-p/po)) (18) 

which reduces to the widely used choice @1j V' — (1 ^ c-^Pi^p)) for po = I. The importance of the free parameter po 
will become apparent later, when discussing the grid-refinement of a given interface. For the moment we confine our 
analysis to the standard case po = 1 and we highlight the role of the two parameters Ai, A2 that, if used properly, allow 
to vary the density separation {pi — Pg)/pg = Ap/ Pg between the two phases and the surface tension a independently 
and in agreement with the continuum interpretation described through (|15p . In figures ([1]) and we show the 
equilibrium profiles obtained with a given Ap/pg at changing the surface tension for both flat and curved surfaces. 
This is done using the two-parameters forcing in such a way to reproduce the same Ai but different A2 in (fTS)) . The 
numerical results for the case of a flat interface are also in good agreement with the theoretical predictions obtained 
from the mechanical stability equation, dyPyy = 0, applied to (|15p . To further check the continuum interpretation 
given through (fTS)) we have carried out Laplace tests for spherical droplets as shown in figure From the Laplace 
law: 



AP = P„, - Pout = 'J/R (19) 

being p„ — Pout the difference between inner and outer pressure in a spherical 2d droplet of radius P, we can estimate 
the surface tension from a lin-lin plot of AP versus 1/P. The numerically estimated surface tension agrees with the 
one predicted by p5|) . pT|) and P?]) . This is the first result presented in this paper. To our knowledge such extension 
of the SC model leading to flexible adjustment of the pressure tensor parameters, i.e. the surface tension and the 
equation of state, has never been considered. This opens the way to describe within a pseudo-potential approach 
more complex physics where the surface tensions needs to be changed independentlyof the equations of state as it is 



the case when, for example, surfactants are added changing the interface properties 421 or when dynamical properties 
must be studied as a function of a as for rising bubbles [43[ . 

In the next subsection we use the extra-freedom given by the tunable reference density po in p8|) to change the 
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numerical resolution of the interface at fixed physics (i.e. fixed surface tension and fixed density ratio). This is an 
important issue, because of the inevitable numerical instabilities which limit the density ratios obtained at a given 
spatial resolution. Indeed, the original SC model ^ is known to be unable to describe density jumps larger than 
Ap/pg ^ 0(10) per grid point. This suggests the possibility of improving the flexibility of the method by spreading 
the same density jump on a larger number of grid nodes. 

Grid refinement and continuum description 

The introduction of po in (fTSl) allows us to refine the interface resolution for fixed density ratio and surface tension. 
If we introduce the shorthand notation: 

V(p) = Vp^^(/5); ^{p) = (1 - exp(-p)) (20) 

where p = p/ po, the pressure tensor ()15|) takes the following expression: 



Pi] = po 



Oid^). (21) 



With reference to the case of a flat interface with dishomogeneities only along the y coordinate, by performing the 
coordinate rescaling y' = Xy in ((2T|) we obtain: 



Ptj = Po 



0{d'^) (22) 



4 ' 2 2 ^' 

where V', A' and 9' means derivatives with respect to the new variable, y' . Let us notice that by choosing 

A2 - A^/A^ (23) 

with A2 constant, the dependency on A disappears from (|22p and the only dependency on po in the above expression 
comes from the overall prefactor. Therefore, in the expression of the mechanical stability condition for a flat interface 
dyPyy = 0, as applied to ([22]), no dependence on po and A is left. In this way, we are able to extract a universal 
profile, V'(p) as a function of y'. This leads to the conclusion that the density ratio Ap/pg is independent of po- 
As for the surface tension, equation ((2T|) in the old variables, yields: 



which, in terms of A'2 and y' , becomes: 



A2P0 
' 2d 



A'2P0 



+ 00 

\dy4'\^dy (24) 

-00 



\dy4\'dy'. (25) 



2Acf 

Since the profile and its integral in the primed variables is universal, from (j25p we see that, by choosing 

A = Po 

the surface tension is also invariant under rescaling of the spatial coordinate. 

It is therefore clear that po in the functional form ([20]) can be used to fine-tune the thickness of the flat interface 
at fixed values of the physical parameters (density ratio and surface tension) provided that we choose A2 = Aj/po- 
In figure we show the equilibrium flat profiles for the case (jTS]) and (pO)) . with Ai = —5.0,^2 = ~5.0/pq and 
different values of po. As one can see the net effect is to change and magnify the interface width with a good agreement 
with the analytical profiles obtained from the continuum description given above. We also carry out (see inset of 
figure ([3])) Laplace tests for the case (fT5|) and ([20]) . with Ai — —5.0, A2 — —5.0/ Pq and three different values of po- The 
macroscopic analysis predicts the same surface tension and indeed this is precisely what the numerical simulations 
show. When moving from large to small po, a refinement of the interface occurs. Thus, fine-tuning of po can be 
regarded as a means of locally magnifying the interface region without changing the macroscopic physics. 
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EQUILIBRIUM DESCRIPTION THROUGH LATTICE BOLTZMANN EQUATIONS 



Up to now, we have mainly investigated the equilibrium properties of interfaces resulting from the addition of a 
pseudo-potential in the classical Lattice Boltzmann formulation. A crucial point is however, to analyze the dynamical 
stability of such results and to understand the effects of the kinematic terms on the equilibrium properties between the 
two phases. For weakly inhomogcncous fluids, this is commonly achieved via the standard Chapman-Enskog expansion 
[Ij Ell using the Knudsen number (molecular mean-free path over smallest macroscopic scale, i.e. the width of the 
interface) as a smallness parameter. However, in the vicinity of a sharp-interface the Knudsen number has no reasons 
to be small, being proportional to density gradients, and the Chapman-Enskog procedure goes under question. Recent 



work in this direction [33|, |3J, |35| has carried out standard Chapman-Enskog analysis with additional forcing terms. 
The proposed analysis leads to a set of different macroscopic dynamic equations. The correctness of the macroscopic 
limit is not analyzed here. Infact, besides detailed analytical control on the behaviour of the hydrodynamic fields 
close to the interface, one may wonder whether numerical implementation of the lattice Boltzmann equation with a 
pseudo-potential provides realistic and stable results over a wide range of density variations and surface tensions. 

Indeed, a disturbing phenomenon, known as spurious currents ^6, 37, Sj, develops systematically in the vicinity 
of interfaces: small circulating currents that are directly proportional to the interface surface tension {i.e. density 
ratio) spoil the physical results of numerical simulations and degrade the numerical stability for high density ratios, 
thus casting serious doubts on the applicability of the method. 



For flat interfaces, the situation is more under control. In fact, all spurious contributions reported near flat interfaces 
are due to an ambiguity in the definition of the fluid momentum. The correct way to measure it, is to take an averaged 
momentum between pre and post-collisional states [s^. 

This cures flat interfaces, but curved interfaces are still affected by the problem and several attempts to justify and 
explain the origin of this phenomenon have been proposed. In (sij . the author proposed an ad-hoc extra counter-term 
to erase spurious currents. Unfortunately, this analysis is limited to flat interfaces and the prescription to erase the 
spurious currents is clearly equivalent to averaging pre and post collisions momentum in the SC model. In [37|, the 
author concluded that the origin of the spurious currents is the incompatibility between the discretization of the 
driving forces for the order parameters and momentum equations. More recently, in [38| . it has been shown that 
spurious currents are due to insufficient isotropy of the discrete forcing operator. In the latter paper, clear numerical 
evidence is brought up, but no detailed analytical explanation is provided. 

Here, besides supporting the numerical findings of [sst . we discuss in details the physical origin of the spurious 
currents. Then, following the symmetry analysis of lattice gas given in (39| . we derive improved isotropic schemes for 
2d and Sd models as well as further possible theoretical improvements. 

The case of flat interface is pretty straightforward. In this case, let us denote again with y the direction of the 
non-homogeneity. We can imagine to have two homogeneous bulk phases p = pg at y = —oo and p = pi at 
y ~ separated by an interface centered aX y = 0. Then, the mass conservation, dtp + dy{puy) = 0, in a 

stationary state {dtp = 0) predicts puy = const, independently of the local density gradients, i.e. independently of 
the Chapman-Enskog expansion Therefore, by imposing a zero net mass- flux at infinity, one readily derives that 
Uy = Q everywhere. 



Let us now analyze the case of a circular drop in 2 dimensions. The new feature is that fluctuations tangential 
to the surface may also appear and their connection with the forcing term plays a key role. Infact, if the forcing is 
perfectly isotropic: 

F,{x) = ^F(r) = e,i?(r), (26) 
r 

where F is a scalar function and is the unit radial vector, one would argue that the velocity field reflects the 
same symmetry, i.e. no spontaneous breaking of rotational invariance should arise. For a stationary state, the 
mass conservation implies: rpuj. = const., being Ur the radial component of the velocity field. The only physical 
acceptable solution is = everywhere. We conclude that if the isotropy of the problem is perfectly carried over 
by the discretization scheme, no spurious currents would develop even for a curved interface. As a consequence, 
the numerically observed currents must be stem from a lack of isotropy at some level with the main contribution to 
anisotropy near the interface due to the pseudo-potential. Indeed, one notices that according to the set of grid points 
and weights entering in the simplest expression of the forcing ([S]), one has a loss of isotropy at a given order in the 
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Taylor expansion. For example, for the simplest case of 2DQ9 one obtains (see Appendix B for details): 

J2 wii\c,\^Wx + ci)ci + + ^V^V^) V. + + + 0{d'), (27) 

where = (1, 0) and = (0, 1) are unit vectors in Cartesian coordinates. If we have axial symmetry of the density 
distribution we must have: "0 = ip(r). Because of V = Bj-dr and = 9^ + r~^dr for an axially symmetric function, 
the first term in the rhs of (|27|1 is isotropic. On the other hand, the 2nd and 3rd terms, that arise only at the fifth 
order, are manifestly anisotropic. We should also notice that in the previous section, we limited our analytical analysis 
to the 4th order expansion, and all the numerical comparison where made by checking that spurious currents arising 
from higher orders were negligible, since we chose a stationary regime with small local gradients in the density field. 
Nevertheless, on the route to higher density ratios, i.e. for cases with high local density gradients, one necessarily 
meets with the problem of anisotropic contributions. In figure (|4|) we show the structure of the spurious currents 
for two cases. As one can see, the currents exhibit typical anisotropics with a quadrupolar modulation, the result 
of anisotropics induced by higher order derivatives in the pseudo-potential expansion (j27[) and they are enhanced 
systematically when the density separation between the two phases is increased. 

To further support the previous statement, we have solved the Laplace equation. Am = with anisotropic boundary 
conditions on a ring, Ur = cos{4:0),ug = for ri < r < r2 (see caption of figure ([5]) for details). The result is a non-zero 
profile in the bulk regions. This is also compared qualitatively with the spurious currents picture from a stationary 
state of a numerical simulation and a good qualitative agreement is observed (see figure (HJ). From these pictures, 
we see that spurious currents, once generated on the interface, spread through the bulk regions, thereby corrupting 
the physical content of numerical simulations. 

Having assessed that spurious currents are triggered by high-order angular harmonics due to lack of sufficient 
isotropy, it is natural to seek new models with a higher degree of isotropy. There are at least two parallel ways to 
remove this problem. Either one improves the support of the underlying lattice structure coupled by the pseudo- 
potential terms, so as to push anisotropy to higher and higher Taylor orders, or one can keep a given degree of 
isotropy of the forcing term and improved grid resolution, so that curved surfaces become more and more refined, 
hence subject to smaller local density gradients. 



Isotropy at a fixed discretization 

The former kind of technical improvement has already been proposed by (38| . Here we support these previous 
findings, and we extend them systematically to higher orders in full details for both 2d and 3d cases (see Appendixes 
C and D). Following [s^, the key idea consists of enlarging only the set of spatial grid points coupled by the pseudo- 
potential tp and choosing the appropriate weights to enforce isotropy up to the desired order. For any practical 
purpose, one writes 

= -gbcl4>ix) wQcil^Mx + ci)c] (28) 

where ci runs over a given set of grid points, changing according to the required order of isotropy (see figure ([6|)) . In 
fact, by applying the Taylor expansion (all details in appendix C) to (|28p . one obtains: 



with 



^^"^ = E^...^ = 5:«;(|c,ncrcr...c- (30) 



(29) 



and (obviously) zero odd tensors: 

E(2n+1) ^ 0. (31) 

The weights can be chosen in such a way to recover isotropy to the desired order (see appendixes C and D) . Clearly, 
more velocities are needed in the implementation of the forcing terms (see figure (|6])). Numerical results (see figure 
^) do confirm a decay of the magnitude of the spurious contributions as the order of isotropy is raised. Although 
the practical implementation of higher-order scheme might not be as straightforward as the standard SC case, it is 
nonetheless reassuring to know that a well-defined procedure to tame spurious currents is available. 
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Refinement at a fixed degree of isotropy 

Since non-isotropic terms in the standard SC forcing scale with fifth-order derivatives, it is plausible to expect 
that these terms can be attenuated also by a refinement of the interface resolution, i.e. by the rescaling procedure 
previously illustrated. In fact, in the standard formulation (|27p the spurious contributions are induced by the terms 

and By that should fade away by a progressive refinement of the grid. 

In figure ([8]) wc show how refining the grid for a fixed surface tension does indeed decrease the amplitude of spurious 
velocities. Using (fT5|) and ([20|) . with the scahng A2 ~ l/Po> macroscopic system stays the same: same surface 
tension and same density ratio. The only difference is a net reduction of the spurious velocity. Let us notice that the 
improvement due to grid-refinement within the extended pseudo-potential (|14p with pressure tensor (jlSp seems more 
effective than the one induced by high-order isotropic forcing in the original SC model. Indeed, comparing figure ([7]) 
and ([8]) one notices that in the latter an almost complete depletion of spurious currents is obtained already with a 
simple factor 2 in the rescaled coordinate. On the other hand, to reach similar level of accuracy in the original SC 
model one needs to improve the isotropy of the forcing up to order 10 or even more. 

The fact that the smoothing of the density profile permits to reduce considerably spurious contributions allows to 
achieve quite large density ratios, up to the order of Ap/ pg ~ 100, as shown in figure where wc plot the maximal 
spurious velocity |t/|maa; normalized to the sound speed as a function of the density ratio. 

Of course, one may also imagine to combine the two proposals, using the extended formulation (|14p with higher 
degrees of isotropy. Whether the numerical effort is worthwhile has to be decided on a casc-by-casc basis. 



CONCLUSIONS 



The SC model is one of the most successful spinoffs of lattice Boltzmann theory. It has nonetheless made the 
object of extensive criticism over the last decade [40j . Part of this criticism is simply misplaced, some other is not. 
In particular, lack of thermodynamic consistency, surface tension ticd-down to the equation of state, and spurious 
currents near sharp interfaces, have spurred doubts on the applicability of the SC method to the simulation of 
realistic multi-phase flows. In this paper we have elucidated the physical reasons behind the above weaknesses, and 
also suggested practical ways around them in the large-scale limit. 

First, we have shown that by enlarging the number of coupling terms in the pseudo-potential expression, one can 
push the method at varying the density ratios and the surface tensions independently and over a wide range of 
parameters. The main limitation in achieving a systematic enhancement of density ratios is due to spurious currents 
in static curved interfaces. This limits both numerical stability in the dynamical evolution and the intimately physical 
correctness even for the static case. 

Second, we have shown how to overcome this problem by developing improved versions of pseudo-potential in- 
teractions. The goal is to reduce anisotropy contributions that are the source of spurious currents. We achieved 
this systematically, either by a refinement of the curved interface, so as to soften the local density gradients, or by 
improving the isotropy of the discretized pseudo-potential. The first method is more effective, leading to a numerical 
reduction of the maximal current up to a factor 10 with only a doubling in the grid resolution. We have shown that 
this stretching of the interface can be achieved by a simple rescaling of the coupling strengths with the reference 
density of the pseudopotential. This permits to achieve an 'adaptive' form of local grid refinement without changing 
the structure of the lattice nodes. 

The present analysis has been carried out for a given choice of the pseudopotential tl'ip)- In principle, the major 
conclusions should carry over to other, possibly more effective, functional forms of tp [30[. 

Besides clarifying the theoretical foundations of the original SC model, it is hoped that the extended version 
presented in this work will help setting the stage for future and more challenging applications of pseudo-potential 
methods to the simulation of complex multiphase flows. 
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APPENDIX A 

In this appendix we discuss the tensorial structures that lead to a vector structure of the form 

J, = ^^jd.Aib = ^^Pdui^p (32) 

where doubled indexes are summed upon. Wc start from the most general expression for a second order, non diagonal 
tensor involving derivatives only in the second order: 



+ d^dui: <5„- + a{^^^P){^J^) + (33) 



Sij 

where a, 6, c, d are meant to be fixed upon consistency with expression (j32p . It is infact verified that upon differenti- 
ation: 

djSij = cdjipdijip + ddiipdjjijj + dipdijjip + adijipdj'ip + adiipdjjip + bdjipdijip + bipdijjip. (34) 
To be consistent with the expression of the forcing wc must impose: 

a+b+c=0 

a + d = (35) 
d + 6= i 

and we end up with three constraints and four constants. This means that there are infinitely many choices of Sij 
satisfying the condition: 

djSij = Ji (36) 

and we need another constraint to close the problem and give unambiguously our tensor. Even if the tensor structure 
is not uniquely determined, when we apply our arguments to the case of a flat interface whose dishomogeneities 
develop along a y coordinate, we notice that the normal component of the above tensor is: 

Syy = l^dyy^ + (« + ^) (dy^)' (37) 
and from the last 2 expressions of (|35|) wc obtain a — b = — i that used in the first one imposes: 

2a + c=-i^a+^ = -i. (38) 
So, even if the tensor is not uniquely determined, its normal component is uniquely given by 

Syy = \'^dyyl}j - '^{dylpf. (39) 

This implies that when using a mechanical stability equation ([TT|) with a fixed boundary condition we are able to 
extract the same profile as a function of y. Then, from the expression (|33p we can also write the equivalent of the 
surface tension considering the mismatch between the normal Syy and tangential components Sxx ■ 

2 



(aidyi'Y + b^dyytl^) = - b) (9^,^) • (40) 

> J —00 

Again, from that the last two expression of (|55l) we get a — b — — ^ . This means that the surface tension is uniquely 
determined. 

APPENDIX B 

In this appendix wc show how to derive non isotropic contributions from discretizations. The forcing term is written 
in the form 

= -GbC^Mx) ^il^il'M^ + ^i)ct- (41) 
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Applying the Taylor expansion, one obtains 



5! 



7! 



(42) 



with 
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(43) 



and (obviously) zero odd tensors: 



£i(2n+l) _ 



(44) 



The even tensors are written as 



where A^^"^ is given by the recursion relation [s 



p(2") _r'(2n)A(2n) 



(45) 



.(4) 



2n 

a(2") = V <5. ■ a(2"-2) f46) 

In our mean field approach, -0 is a function of the density. If the density distribution is axially symmetric, ip is also 
axially symmetric, ip = ^(^)- Then, the force Fi should be written as 

F,ix) = ^F{r), (47) 
r 

where F is a scalar function. It should be noted that the isotropy for all £;(2n) jg essential in order to satisfy the 
relation l|47p . Now, we will show that the truncated isotropy induces the anisotropic force, which triggers the spurious 
currents, even when the density distribution is axially symmetric. Let us consider the standard case of 2DQ9. As 
already noticed in the text, this is a 4th-ordcr approximation in the isotropy and the weights are given by 

«;(l) = l/3, u;(2) = l/12. 



w{n) = for n > 3. (48) 

This approximation means that all the tensors up to the 4th-ordcr (£j'^' and iJ'^'*)) arc isotropic but the higher order 
ones (n > 6) are not. Using standard Taylor expansion for lattice Boltzmann populations one obtains after some 
lengthly algebra: 

J2 HWil'm^ + ci)c[ + + d^)^ + (^^dt + ^dld^ + ^d^^ ^ + 

' / \ ^^^^ 

Using a nabla operator, is rewritten as 

^ w{\ci\^mx + ci)ci ^v(l + \v^ + ^V^V^) + + e,^ + 0{d'), (50) 

where = (1,0) and By = (0, 1) are unit vectors in Cartesian coordinates. Next we assume axial symmetry of the 
density distribution, i.e., if} — tp{r). Because of V = Bj-dr and = + r~^dr for an axially symmetric function. 
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the 1st term in the r.h.s of ((50|) is isotropic. On the other hand, the 2nd and 3rd terms, related to the anisotropic 
tensor E'^^'^ in (|42p . are not. Now the force Fi is decomposed into the isotropic and anisotropic parts, i.e., 

F,{x) = ^F{r) + F^{x). (51) 

Within the 0{d^) approximation, F(r) and F({x) arc respectively given by 

Idd lddidd\ 

^ 180 ' ^ 180 ' ^ ' 

The anisotropic force F/ due to the anisotropy of E''^^ is responsible for the spurious currents. Higher orders can be 
computed similarly (for the interested reader please contact the authors). 



APPENDIX C 



Here we detail the exact procedures leading to higher order isotropic terms in the forcing contribution for a regular 
lattice in 2d. The velocity phase space and forcing weights for isotropic terms up to 16th order are explicitly given 
in figure To treat correctly isotropy from a lattice set of velocity vectors [c\,i = x,y) the starting point is the 2 
point tensor on the lattice which is assumed normalized to unity 

Y,<c^MWi?)=5,3- (54) 
I 

Considering the regular structure of the lattice and the consequent symmetry of c\ with respect to i = x and i = y, 
one can write (|54p in the simplified form 

^(cDMl^n-l. (55) 

The fourth-order isotropy is imposed by 

J2 c^c^cXIqI^) = C(4) (%4s + 6,kS,s + , (56) 

where C*-^-* is a constant. Since one obtains 

^(cnM|c|?) = 3CW, ^(cn^(cF)M|cd^)=C(^\ 

a condition to satisfy (|56|) is written as 



Eiicf)M\ci\') 

j:iicmcf)M\ci\' 



3. (57) 



In terms of lattice vector this can be achieved with the standard 2DQ9 model with weights w(l) = 1 /3 and w(2) = 1/12 

and the corresponding lattice velocities: 

w{l): 

(Ci C2 C3 C4) = g ^ -1 ) ' ^^^"l 

w{2): 

(C5 Cg C7 Cg) = ^ I 1 -1 -1 ) • 
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More general conditions can then be obtained for higher order tensors. For example, the 6th. 8th, 10th and 12th-order 
isotropics are given by 

,3 , 10 

/ 

^ ^ U7(|c/| ) — ^ ('^iii2 '^i.sM ^^^5*6 ^^7^8 ^^9^10 '^^11^12 ^~ ■') • 

And the mixed contributions can be constructed as well: 

^(cf )'"(cf)'"i«(|Q|2) = c(2"+2")(2n - l)!!(2m - 1)!!, (61) 

where (2n — 1)!! = (2?! — 1) x (2?! — 3) x ... x 1. Then, to achieve isotropy at higher orders one should introduce some 
requirements on the tensors. Just to give an example, for the isotropy up to the 6th order one should require that: 



i:i{c^ifH\ci\^) 

E,(cr)4(cf)2«;(|c,P) 



(62) 



5. (63) 



these translate to the matrix relation 



2 4 8 \ (w{l)\ (1 

1 -4 16 w{2) = I (64) 
1 -8 64 / \ w(4) ) \ 

that can be satisfied using 12 velocities with weights w(l) = 4/15, w(2) = 1/10 and w(4) = 1/120 
w{l): 

(Ci C2 C3 C4) = ( J ^ -1 ' ' ^^^^ 



w{2): 



w(4): 



(C5 C6 C7 Cg) = ( J 1 -1 -1 j ' 



/ ,,20-20, 

(cg Cio Cn C12) = ( Q 2 -2 ' ■ -* 



Higher order calculations are lengthly and not reported here. The set of vectors can be extracted from figure ^ while 
the weights can be found in table 1. 



APPENDIX D 



The same calculations are then arranged in id. For each w{n) (reported in table 2), the corresponding velocity 
vectors c; are shown below: 



w{l): 



w{2): 



w{3): 



w(4): 



w{5): 



w{6): 



w{8): 



W22l(9) 



W3oo(9) 



/ 1 -1 \ 

(Ci C2 Cg C4 Cg Cg) = 1 -1 , 

\ 1 -1 / 



(C7 Cg Cg Cio Cii Ci2 Ci3 C14 C15 Cig C17 Cig) 

1 
1 





1 


-1 


-1 


1 


1 


-1 


-1 











1 


1 


-1 














1 1 


-1 


-1 











1 


-1 


1 


-1 


1 -1 


1 


-1 



/ 1 1 1 1-1-1 

(Ci9 Cio C21 C22 C23 C24 C2,5 C26) =1 1-1-1 1 1 

\ 1 -1 1-1 1 1 



/2-20 00 o\ 

(C27 C28 C29 C30 C31 C32) = 2 -2 

\0 00 02 -2 / 



(C33 C34 C35 C36 C37 C38 C39 C40 C41 C42 C43 C44 

C45 C46 C47 C48 C49 C50 C51 C52 C53 C54 C,55 Cgg) 





2 




2 




2 2 


2 




2 


-2 1 


-1 




1 


(; 


-1 




1 




1 










2 


2 




2 














1 


-1 




1 


-1 
























1 - 


-1 


1 




1 










2 


2 - 


-2 




2 













1 - 


-1 


1 




1 - 


-1 


1 




1 


2 


2 - 


-2 




-2 2 


2 - 


-2 





(Csi .. C92) 



2 .. 


. -1 


1 .. 


. -1 


1 .. 


. -2 


2 .. 


. 


2 .. 


. -2 


.. 


. -2 



(C93 .. Cue) = 



(Cii7 .. C122) 
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w(10): 



(ci23 .. Cms) 



w{ll): 



(Ci47 .. Ci7o) 



(77) 



(78) 
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FIG. 1; Smoothing of interface properties. We show the static profiles for a flat interface with an improved SC model (|14|l and 
(I7|). The parameters are chosen to produce a macroscopic pressure tensor (|15|l keeping fixed Ai — —5.0 (same density ratio) 
and at varying A2: A2 = —5.0 (□) and A2 = —30.0 (o). The results of numerical simulations are compared with the analytical 
estimates (solid lines) obtained solving the mechanical stability equation pip applied to (|15p . Notice the smoothing in the 
interface for a fixed density ratio due to a change in the surface tension. In all numerical cases the lattice Boltzmann equation 
U]) has been integrated in time in a fully periodic domain x Ly = lOOAa; x lOOAa; with a fiat strip of liquid in an otherwise 
gaseous domain. The relaxation time is r = 0.7 At. 
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FIG. 2: Laplace test for the same density ratio obtained from the numerical simulation of the improved SC model (|14p and ((Tj. 
The parameters are chosen to reproduce a macroscopic pressure tensor (|15|) with Ai = —5.0 and different A2: A2 = —5.0 (□), 
A2 = —15.0 (o), A2 = —30.0 (A). For each case we plot the pressure drop AP as a function of the inverse radius of the static 
drop. Moreover we compare the results with the theoretical predictions (solid lines) given by the continuum analysis that leads 
to a{A2 — —5) = 0.0398, a{A2 = —15) = 0.0716 and a{A2 — —30) — 0.100 in lattice units. Numerical details are the same of 
figure U]) with the only difference that now we use a drop in the middle of the domain. 
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FIG. 3: Smoothing of interface properties. We show the static profiles for a flat interface with an improved SC model (|14|l and 
(|20|) . The parameters are chosen such as to produce a macroscopic pressure tensor (|2H) with Ai — —5.0 and A2 — — 5.0/po for 
different values of po: po = 1.0 (□), po = 0.5 (o) . Results are also compared with the analytical estimate resulting from our 
continuum interpretation (solid line). The two profiles have been plotted by rescaling the lattice grid by a factor 1/po- Notice 
the increased interface grid resolution obtained at fixed density ratio. Inset: the results of the Laplace tests made on spherical 
droplets with the same density ratio and varying grid resolution is also plotted. Again, we get the same surface tension with 
different interface resolutions, po = 1.0 (□), po ~ 0.7 (o), po — 0.5 (A). In all numerical cases the lattice Boltzmann equation 
has been integrated in time in a fully periodic domain x Ly — lOOAa; x lOOAa: with a flat strip of liquid in an otherwise 
gaseous domain. The relaxation time is r = 0.7 At. 






FIG. 4: Spurious contributions around a static 2d drop for different density ratios with the standard SC model ([6]) and ([7|. 
Left: we show the local Mach number defined as -f /cs with Cs the lattice sound speed. A drop of radius I5A2; and 
density ratio Ap/pg ~ 35, obtained with Qb = —6.0 in (|6]), is considered. Right: the same as in the case of left figure with a 
higher density ratio Ap/pg ~ 60 obtained with Qt ~ —7.0 in ([6]). Note in both plots the angular dependency due to lack of 
perfect radial symmetry in the forcing terms. In all numerical cases the lattice Boltzmann equation ([1} has been integrated in 
time in a fully periodic 2d domain Lx x Ly = lOOAa; x lOOAx with the drop put in the middle of the system. The relaxation 
time is r = 1.0 At. 
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FIG. 5: Spurious contributions in Lattice Boltzmann and their continuum interpretation. Left: we report the structure of the 
velocity field around a static drop of radius 15Ax and density ratio Ap/pg ~ 50. The data are the same of right panel of 
figure ([4]). Right: for a qualitative comparison we have solved the Laplace equation together with the continuity equations 
with a matching condition in a radial ring of a given width. The velocity field is obtained in an iterative way by first solving 
the Poisson equation = V ■ u and then renewing the velocity filed as u — > u — Vfj). For both u and (f) periodic boundary 
conditions are imposed in the horizontal and vertical directions. On the ring of range 17.5Aa; < r < 20.5Aa; the matching 
condition: Ur = cosAO, Me = is imposed being Ur and ug the radial and azimuthal components of the velocity field respectively. 




FIG. 6: The grid points identifying the set of velocity fields for a 2d case. With reference to the weights reported in table 
1, different degrees of isotropy can be achieved: 4th order (up to w{2)), 6th order (up to w(4)), 8th order (up to ?i'(8)), 10th 
order (up to w(10)), 12th order (up to w(17)), 14th order (up to U)(25)) and 16th order (up to w[2>2)). 
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FIG. 7; Reduction of spurious currents with higher order isotropic tensors in the forcing terms. The spurious contributions for 
a static drop are analyzed for a density ratio Ap/pj ~ 50 obtained for a standard SC model ((Gjl and ((7|). The parameter chosen 
is Gb = —7.0. We show the vertical velocity, Uy, normalized with the lattice speed of sound, Ca, at x = Lx/2 and as a function 
of y/ Ax. The different plots correspond to different degrees of isotropy in the forcing term ([S}: 4th order (□), 6th order (o), 
8th order (A) and 10th order (x). Notice the reduction of the spurious contributions in the proximity of the droplet surface, 
y/Ax ~ 30, 70. Data are obtained from Lattice Boltzmann equation ((T} in a fully periodic domain x Ly — lOOAx x lOOAi 
and T — l.OAt in lattice units. 
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FIG. 8: reduction of the spurious currents with grid refinement at changing po. The spurious contributions for a static drop 
are analyzed for a density ratio Ap/pg ~ 50 obtained for the improved SC model (|14|l and (|20p . The parameters are chosen in 
such a way to reproduce (|15p with Ai = —7.0. We show the vertical velocity, Uy, normalized with the lattice speed of sound, 
Cs, at a; = Lx/2 ad as a function of y/Ax. The different plots correspond to different degrees of refinement obtained with 
A2 = — 7.0/po in PSI) : po = 1.0 (□), po — 0.75 (o) and po = 0.5 (A). Details of the numerical simulations are the same of the 
previous figure. 
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FIG. 9: Reduction of the spurious currents with refinement in the forcing terms. The spurious contributions for a static drop 
are analyzed for various density ratios obtained in the improved SC model (|14|l and (|20|) . The parameters are chosen in such 
a way as to reproduce (|15[) with different Ai, thus changing the density ratio Ap/ pg. We show the maximum velocity due to 
spurious contributions, normalized with the lattice speed of sound (cs) as a function of Ap/pg. The 2 different plots correspond 
to different degree of refinement obtained with A2 = 7.0/po in (EI} a-nd different values of po: po = 1.0 (□), po = 0.5 (o). 
Notice the net reduction of the spurious contributions at fixed macroscopic physics (same density ratio and surface tension) 
obtained through the rescaling A2 ~ 1/po- The numerical details are the same described in figure ((Sjl. 
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Tabic 1. Weights up to the 16th-order approximation for the case of 2d models. Notice that the weights for 
velocities with |c/p = 25 needs to be chosen differently according to the directions in the x — y plane. The notation 



') stands for the velocity lattice vectors (±a;i,±a;2) and (±a;2,±a;i). 
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^(10) 
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w(l) 
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4 
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Table 2 Weights up to the lOth-order approximation for 3d models. Notice that the weights for velocities with 
|qP = 9 needs to be chosen differently according to the directions in the x ~ y — z space. The notation 'Wx,y.z{\ci\'^) 
stands for the velocity lattice vectors (±a;i, ±X2, ±3:3) plus permutation. 
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